package org.zhongyin.math;

/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */


/**
 *
 * @author yguo4
 */
import java.io.*;
public class Regression
{

    /**
     * @param args the command line arguments
     */
    public static void main(String[] args)
    {
        di nput=new di();
        mr24 mrs=new mr24();
//        mf ore=new mf();
        double[][] a=nput.inp("c:\\temp\\Test.txt");
//        System.out.println(a[0][0]);
        double[][] d=mrs.reg(a);
        int row=d[0].length;
        for (int i=0;i<row;i++)
        {
            //System.out.print("--"+i+"--"+"\n");
            for (int j=0;j<48;j++)
            {
                System.out.printf("%11.6f"+"\t",d[j][i]);
            }  
            System.out.print("\n");
        }


    }
}





//==============================================================================
//define the class for data loading
class di
{
//------------------------------------------------------------------------------      
    //define the # of replicates from original input
    private int rowinp;
    //define the # of variables from original input
    private int colinp;
    //char reader for data reading
    private char chinp;
    //temporary string saver for converting data from string to double
    private String tempinp;
    //define the output of data reading
    private double[][] xinp;
    
    //data loading from local file
    double[][] inp(String str)
    //input: location of local file
            
    //output: data from local file,
    //Column 1: month; Column 2: day of week; Column 3: hour of day; Column 4: cooling load at time t;
    //Column 5: daily maximum temperature; Column 6: daily minimum temperature; Column 7: temperature at time t;
    //Column 8: relative humidity at time t; Column 9: data type indicator, 0 for regression, 1 for forecasting
            
    {
        FileInputStream fin;
        DataInputStream din;
        //define the initial value of # of replicates
        rowinp=1;
        //define the # of column
        colinp=9;
        //count the # of replicates in the data of local file
        try
        {
            fin=new FileInputStream(str);
            din=new DataInputStream(fin);
            while((chinp=din.readChar())!='\0')
            {
                if(chinp=='\n')
                //a line is one replicate
                {
                    rowinp=rowinp+1;
                }
            }
            din.close();
            fin.close();
        }
        catch(FileNotFoundException e)
        {
            System.out.println("File Not Found!");
        }
        catch(IOException e)
        {
                //^_^//
        }
        //define the matrix of data to be loaded
        xinp=new double[rowinp][colinp];
        try
        {
            fin=new FileInputStream(str);
            din=new DataInputStream(fin);
            //initialize string buffer to be empty
            tempinp="";
            //define the index of row and column for the data matrix
            int iinp=0;
            int jinp=0;
            //read char until data ends
            while((chinp=din.readChar())!='\0')
            {
                //throw away invalid char
                if(chinp!='﻿')
                {   
                    //load single number and assign to the data matrix
                    //tab, space, new line are the separator character
                    if(((chinp=='\t' | chinp==' ') | chinp=='\n') & !tempinp.equals(""))
                    {
                        //convert string to double
                        xinp[iinp][jinp]=Double.parseDouble(tempinp);
                        //move index
                        jinp=jinp+1;
                        if(jinp==colinp)
                        {
                            jinp=0;
                            iinp=iinp+1;
                        }
                        //initialize string buffer to be empty for loading next value
                        tempinp="";
                    }
                    //remove repeating separator charactors
                    else if(((chinp=='\t' | chinp==' ') | chinp=='\n') & tempinp.equals(""))
                    {
                      //WoW//
                    }
                    else
                    {
                        //read valid char and add to the end of string
                        tempinp=tempinp+chinp;
                    }
                }
            }
            din.close();
            fin.close();
        }
        catch(FileNotFoundException e)
        {
            System.out.println("File Not Found!");
        }
        catch(IOException e)
        {
                //ToT//
        }
        return xinp;
    }
}





//==============================================================================
//define the class for multiple regression of 24-hour model
class mr24 
{
    //--------------------------------------------------------------------------
    //define the # of replicates
    private int rowinppre;
    //define the # of different types of value
    private int colinppre;
    //define the index of # of replicates for each hour
    private int[] ninppre;
    //define the value of hour hour 1-24 are for weekday, hour 25-48 are for weekend
    private int hourinppre;
    //define the matrix for processed input data
    private double[][][] xinppre;
    
    //input data preprocess
    double[][][] reginppre24(double[][] array)
    //input: array, original data;
    //Column 1: month; Column 2: day of week; Column 3: hour of day; Column 4: cooling load at time t;
    //Column 5: daily maximum temperature; Column 6: daily minimum temperature; Column 7: temperature at time t;
    //Column 8: relative humidity at time t; Column 9: data type indicator, 0 for regression, 1 for forecasting
            
    //output: X & Y when index=0 for regression, X & Y=0 when index=1 for forcasting
    //Column 1: cooling load at time t-24; Column 2: cooling load at time t-168; Column 3: daily maximum temperature;
    //Column 4: daily minimum temperature; Column 5: temperature at time t; Column 6: temperature at time t-24;
    //Column 7: temperature at time t-168; Column 8: relative humidity at time t;    
    //Column 9: relative humidity at time t-24; Colmun 10: relative humidity at time t-168;
    //Column 10: cooling load at time t
            
    {
        //obtain or define the demension of output matrix
        rowinppre=array.length;
        colinppre=11;
        //define the matrix of data after data loading and input data preprocess
        xinppre=new double[48][rowinppre/24][colinppre];
        //define the index for the data of each hour
        ninppre=new int[48];
        
        //data sorting for regression
        //start sorting all the data, first 168 data points need to be reserved for time lag reference
        for (int i=168;i<rowinppre;i++)
        {                //check the index for regression
            if (array[i][8]==0)
            {
                //data for regression must between Apr and Oct
                if(array[i][0]>3 & array[i][0]<11)
                {
                    //insure the last two data points will not be reached when performing data checking
                    if(i<rowinppre-2)
                    {
                        //monday
                        if(array[i][1]==1)
                        {
                            //obnormal reading type one: obnormal peak or valley value
                            if((array[i-1][3]<2*array[i][3] & array[i+1][3]<2*array[i][3]) & (array[i-1][3]>0.5*array[i][3] & array[i+1][3]>0.5*array[i][3]))
                            {
                                if((array[i-1-72][3]<2*array[i-72][3] & array[i+1-72][3]<2*array[i-72][3]) & (array[i-1-72][3]>0.5*array[i-72][3] & array[i+1-72][3]>0.5*array[i-72][3]))
                                {
                                    if((array[i-1-168][3]<2*array[i-168][3] & array[i+1-168][3]<2*array[i-168][3]) & (array[i-1-168][3]>0.5*array[i-168][3] & array[i+1-168][3]>0.5*array[i-168][3]))
                                    {
                                        //abnormal reading type two: null value
                                        if((array[i][3]!=0 & array[i-72][3]!=0) & array[i-168][3]!=0)
                                        {
                                            //obnormal reading type three: continuous same reading
                                            if((array[i][3]!=array[i-1][3] & array[i][3]!=array[i-2][3]) & ((array[i][3]!=array[i-1][3] & array[i][3]!=array[i+1][3]) &(array[i][3]!=array[i+2][3] & array[i][3]!=array[i+1][3])))
                                            {
                                                if((array[i-72][3]!=array[i-1-72][3] & array[i-72][3]!=array[i-2-72][3]) & ((array[i-72][3]!=array[i-1-72][3] & array[i-72][3]!=array[i+1-72][3]) &(array[i-72][3]!=array[i+2-72][3] & array[i-72][3]!=array[i+1-72][3])))
                                                {
                                                    if((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i-2-168][3]) & ((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i+1-168][3]) &(array[i-168][3]!=array[i+2-168][3] & array[i-168][3]!=array[i+1-168][3])))
                                                    {
                                                        //obtain the value of hour
                                                        hourinppre=(int)array[i][2];
                                                        //time lag reference, for Monday one day back is last Friday
                                                        //cooling load at t-24 (one day)
                                                        xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-72][3];
                                                        //cooling load at t-168 (one week)
                                                        xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                                        //maximum temperature
                                                        xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                                        //minimum temperature
                                                        xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                                        //temperature at time t
                                                        xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                                        //temperature at time t-24 (one day)
                                                        xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-72][6];
                                                        //temperature at time t-168 (one week)
                                                        xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                                        //relative humidity at time t
                                                        xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                                        //relative humidity at time t-24 (one day)
                                                        xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-72][7];
                                                        //relative humidity at time t-168 (one week)
                                                        xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                                        //cooling load at time t (dependent variable)
                                                        xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                                        //index increase by 1
                                                        ninppre[hourinppre]=ninppre[hourinppre]+1;
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                        //Tuesday to Friday
                        else if(array[i][1]<6)
                        {                            
                            if((array[i-1][3]<2*array[i][3] & array[i+1][3]<2*array[i][3]) & (array[i-1][3]>0.5*array[i][3] & array[i+1][3]>0.5*array[i][3]))
                            {
                                if((array[i-1-24][3]<2*array[i-24][3] & array[i+1-24][3]<2*array[i-24][3]) & (array[i-1-24][3]>0.5*array[i-24][3] & array[i+1-24][3]>0.5*array[i-24][3]))
                                {
                                    if((array[i-1-168][3]<2*array[i-168][3] & array[i+1-168][3]<2*array[i-168][3]) & (array[i-1-168][3]>0.5*array[i-168][3] & array[i+1-168][3]>0.5*array[i-168][3]))
                                    {
                                        if((array[i][3]!=0 & array[i-24][3]!=0) & array[i-168][3]!=0)
                                        {
                                            if((array[i][3]!=array[i-1][3] & array[i][3]!=array[i-2][3]) & ((array[i][3]!=array[i-1][3] & array[i][3]!=array[i+1][3]) &(array[i][3]!=array[i+2][3] & array[i][3]!=array[i+1][3])))
                                            {
                                                if((array[i-24][3]!=array[i-1-24][3] & array[i-24][3]!=array[i-2-24][3]) & ((array[i-24][3]!=array[i-1-24][3] & array[i-24][3]!=array[i+1-24][3]) &(array[i-24][3]!=array[i+2-24][3] & array[i-24][3]!=array[i+1-24][3])))
                                                {
                                                    if((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i-2-168][3]) & ((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i+1-168][3]) &(array[i-168][3]!=array[i+2-168][3] & array[i-168][3]!=array[i+1-168][3])))
                                                    {
                                                        hourinppre=(int)array[i][2];
                                                        //time lag reference, for Tuesday to Friday is yesterday
                                                        xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-24][3];
                                                        xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                                        xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                                        xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                                        xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-24][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-24][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                                        ninppre[hourinppre]=ninppre[hourinppre]+1;    
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                        else if(array[i][1]==6)  
                        {
                            if((array[i-1][3]<2*array[i][3] & array[i+1][3]<2*array[i][3]) & (array[i-1][3]>0.5*array[i][3] & array[i+1][3]>0.5*array[i][3]))
                            {
                                if((array[i-1-144][3]<2*array[i-144][3] & array[i+1-144][3]<2*array[i-144][3]) & (array[i-1-144][3]>0.5*array[i-144][3] & array[i+1-144][3]>0.5*array[i-144][3]))
                                {
                                    if((array[i-1-168][3]<2*array[i-168][3] & array[i+1-168][3]<2*array[i-168][3]) & (array[i-1-168][3]>0.5*array[i-168][3] & array[i+1-168][3]>0.5*array[i-168][3]))
                                    {
                                        if((array[i][3]!=0 & array[i-144][3]!=0) & array[i-168][3]!=0)
                                        {
                                            if((array[i][3]!=array[i-1][3] & array[i][3]!=array[i-2][3]) & ((array[i][3]!=array[i-1][3] & array[i][3]!=array[i+1][3]) &(array[i][3]!=array[i+2][3] & array[i][3]!=array[i+1][3])))
                                            {
                                                if((array[i-144][3]!=array[i-1-144][3] & array[i-144][3]!=array[i-2-144][3]) & ((array[i-144][3]!=array[i-1-144][3] & array[i-144][3]!=array[i+1-144][3]) &(array[i-144][3]!=array[i+2-144][3] & array[i-144][3]!=array[i+1-144][3])))
                                                {
                                                    if((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i-2-168][3]) & ((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i+1-168][3]) &(array[i-168][3]!=array[i+2-168][3] & array[i-168][3]!=array[i+1-168][3])))
                                                    {
                                                        hourinppre=(int)array[i][2]+24;
                                                        //time lag reference, for Saturday is last Sunday
                                                        xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-144][3];
                                                        xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                                        xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                                        xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                                        xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-144][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-144][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                                        ninppre[hourinppre]=ninppre[hourinppre]+1;
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                        else
                        {
                            if((array[i-1][3]<2*array[i][3] & array[i+1][3]<2*array[i][3]) & (array[i-1][3]>0.5*array[i][3] & array[i+1][3]>0.5*array[i][3]))
                            {
                                if((array[i-1-24][3]<2*array[i-24][3] & array[i+1-24][3]<2*array[i-24][3]) & (array[i-1-24][3]>0.5*array[i-24][3] & array[i+1-24][3]>0.5*array[i-24][3]))
                                {
                                    if((array[i-1-168][3]<2*array[i-168][3] & array[i+1-168][3]<2*array[i-168][3]) & (array[i-1-168][3]>0.5*array[i-168][3] & array[i+1-168][3]>0.5*array[i-168][3]))
                                    {
                                        if((array[i][3]!=0 & array[i-24][3]!=0) & array[i-168][3]!=0)
                                        {
                                            if((array[i][3]!=array[i-1][3] & array[i][3]!=array[i-2][3]) & ((array[i][3]!=array[i-1][3] & array[i][3]!=array[i+1][3]) &(array[i][3]!=array[i+2][3] & array[i][3]!=array[i+1][3])))
                                            {
                                                if((array[i-24][3]!=array[i-1-24][3] & array[i-24][3]!=array[i-2-24][3]) & ((array[i-24][3]!=array[i-1-24][3] & array[i-24][3]!=array[i+1-24][3]) &(array[i-24][3]!=array[i+2-24][3] & array[i-24][3]!=array[i+1-24][3])))
                                                {
                                                    if((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i-2-168][3]) & ((array[i-168][3]!=array[i-1-168][3] & array[i-168][3]!=array[i+1-168][3]) &(array[i-168][3]!=array[i+2-168][3] & array[i-168][3]!=array[i+1-168][3])))
                                                    {
                                                        hourinppre=(int)array[i][2]+24;
                                                        //time lag reference, for Sunday is yesterday
                                                        xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-24][3];
                                                        xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                                        xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                                        xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                                        xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-24][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                                        xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-24][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                                        xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                                        ninppre[hourinppre]=ninppre[hourinppre]+1;
                                                    }
                                                }
                                            }
                                        }
                                    }
                                }
                            }
                        }
                    }
                    //last two data points will only check type two obnormal reading 
                    else
                    {
                        if (array[i][1]==1)
                        {
                            if ((array[i][3]!=0 & array[i-72][3]!=0) & array[i-168][3]!=0)
                            {
                                hourinppre=(int)array[i][2];
                                xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-72][3];
                                xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-72][6];
                                xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-72][7];
                                xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                ninppre[hourinppre]=ninppre[hourinppre]+1;
                            }
                        }
                        else if(array[i][1]<6)
                        {
                            if ((array[i][3]!=0 & array[i-24][3]!=0) & array[i-168][3]!=0)
                            {
                                hourinppre=(int)array[i][2];
                                xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-24][3];
                                xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-24][6];
                                xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-24][7];
                                xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                ninppre[hourinppre]=ninppre[hourinppre]+1;
                            }
                        }
                        else if(array[i][1]==6)
                        {
                            if((array[i][3]!=0 & array[i-144][3]!=0) & array[i-168][3]!=0)
                            {
                                hourinppre=(int)array[i][2]+24;
                                xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-144][3];
                                xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-144][6];
                                xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-144][7];
                                xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                ninppre[hourinppre]=ninppre[hourinppre]+1;
                            }   
                        }
                        else
                        {
                            if((array[i][3]!=0 & array[i-24][3]!=0) & array[i-168][3]!=0)
                            {
                                hourinppre=(int)array[i][2]+24;
                                xinppre[hourinppre][ninppre[hourinppre]][0]=array[i-24][3];
                                xinppre[hourinppre][ninppre[hourinppre]][1]=array[i-168][3];
                                xinppre[hourinppre][ninppre[hourinppre]][2]=array[i][4];
                                xinppre[hourinppre][ninppre[hourinppre]][3]=array[i][5];
                                xinppre[hourinppre][ninppre[hourinppre]][4]=array[i][6];
                                xinppre[hourinppre][ninppre[hourinppre]][5]=array[i-24][6];
                                xinppre[hourinppre][ninppre[hourinppre]][6]=array[i-168][6];
                                xinppre[hourinppre][ninppre[hourinppre]][7]=array[i][7];
                                xinppre[hourinppre][ninppre[hourinppre]][8]=array[i-24][7];
                                xinppre[hourinppre][ninppre[hourinppre]][9]=array[i-168][7];
                                xinppre[hourinppre][ninppre[hourinppre]][10]=array[i][3];
                                ninppre[hourinppre]=ninppre[hourinppre]+1;
                            }
                        }
                    }    
                }
            }
        }
        return xinppre;
    }
//------------------------------------------------------------------------------      
    //define the matrix for the preprocess of independent variables
    private double[][] xpre;
    //define the # of input, # of independent variables, # of replicates
    private int rowpre;
    private int colpre;
    private int npre;
    
    //independent variable preprocess
    double[][] regpre24(double[][][] array,int index)
    //input: array, loaded data from local file; index for the day of week
    //Column 1: cooling load at time t-24; Column 2: cooling load at time t-168; Column 3: daily maximum temperature;
    //Column 4: daily minimum temperature; Column 5: temperature at time t; Column 6: temperature at time t-24;
    //Column 7: temperature at time t-168; Column 8: relative humidity at time t;    
    //Column 9: relative humidity at time t-24; Colmun 10: relative humidity at time t-168;
    //Column 10: cooling load at time t        
            
    //output: processed variables
    //Column 1: constant (1); Column 2: cooling load at time t-24; Column 3: cooling load at time t-168;
    //Column 4: daily maximum temperature; Column 5: daily minimum temperature; Column 6: temperature at time t;
    //Column 7: cooling load at time t-24; Column 8: cooling load at time t-168; Column 9: relative humidity at time t;
    //Column 10: relative humidity at time t-24; Column 11: relative humidity at time t-168;
    //Column 12: temperature square at time t; Column 13: temperature square at time t-24;
    //Column 14: temperature square at time t-168; Column 15: product of temperature and relative humidity at time t;
    //Column 16: product of temperature and relative humidity at time t-24;
    //Column 17: product of temperature and relative humidity at time t-168;
    //Column 18: cooling load at time t (dependent variable)        
            
    {
        //get the row (number of replicates) of original data
        rowpre=array[index].length;
        //get the column (number of orignal dependent variables) of original data
        colpre=array[index][0].length; //=11
        for (int i=0;i<rowpre;i++)
        {
            if(array[index][i][0]==0 & array[index][i][1]==0 & array[index][i][2]==0)
            {
                rowpre=i;
                break;
            }
        }
        //the number of variables (n): constant(Col 1)+independent(Col 2-17)+dependent(Col 18)
        npre=1+colpre+3+3; //=18
        //define the matrix of output
        xpre=new double[rowpre][npre];
        //data preprocess
        for (int i=0;i<rowpre;i++)
        {
            for (int j=0;j<npre;j++)
            {
                //setup constant term
                if (j==0)
                {
                    xpre[i][j]=1;
                }
                //copying the original independent variables
                else if (j<colpre)
                {
                    xpre[i][j]=array[index][i][j-1]; 
                }
                //generate the higher order term from original independent variables
                else if (j<colpre+3)
                {
                    xpre[i][j]=array[index][i][j-7]*array[index][i][j-7];
                }
                //generate the interaction term from original independent variables
                else if (j<colpre+6)
                {
                    xpre[i][j]=array[index][i][j-10]*array[index][i][j-7];
                }
                else
                {
                    xpre[i][j]=array[index][i][j-7];
                }
            }
        }
        return xpre;
    }
//------------------------------------------------------------------------------    
    //define the matrix for initial regression
    private double[][] xreginit;
    private double[][] yreginit;
    private double[][] betareginit;
    //define the # of independent variable for initial regression
    private int colreginit;
    //define the # of replicates
    private int rowreginit;
    //define the intermediate variable for initial regression calculation
    private double[][] xtinit;
    private double[][] xtxinit;
    private double[][] xtxinvinit;
    private double[][] xtxinvxtinit;
    
    //multiple regression
    double[][] reginit(double[][] array1,double[][] array2)
    //input: array1, independent variable; array2, dependent variable
    //Column 1: constant (1); Column 2: cooling load at time t-24; Column 3: cooling load at time t-168;
    //Column 4: daily maximum temperature; Column 5: daily minimum temperature; Column 6: temperature at time t;
    //Column 7: cooling load at time t-24; Column 8: cooling load at time t-168; Column 9: relative humidity at time t;
    //Column 10: relative humidity at time t-24; Column 11: relative humidity at time t-168;
    //Column 12: temperature square at time t; Column 13: temperature square at time t-24;
    //Column 14: temperature square at time t-168; Column 15: product of temperature and relative humidity at time t;
    //Column 16: product of temperature and relative humidity at time t-24;
    //Column 17: product of temperature and relative humidity at time t-168;
            
    //output: beta
            
    {
        //function call for matrix operation
        mxop erareginit=new mxop();
        //obtain the demension of input matrix
        rowreginit=array1.length;
        colreginit=array1[0].length; //=17
        //define new matrix for temporarily saving input data
        xreginit=new double[rowreginit][colreginit];
        yreginit=new double[rowreginit][1];
        //copy input data from input matrix
        for (int i=0;i<rowreginit;i++)
        {
            //copy the independent variables
            System.arraycopy(array1[i], 0, xreginit[i], 0, colreginit);
            //copy the dependent variable
            System.arraycopy(array2[i], 0, yreginit[i], 0, 1);
        }
        //define beta
        betareginit=new double[colreginit][1];
        //calculate beta=(xT*x)^(-1)*xT*y
        //xT
        xtinit=erareginit.trans(xreginit);
        //xT*x
        xtxinit=erareginit.mult(xtinit,xreginit);
        //(xT*x)^(-1)
        xtxinvinit=erareginit.inv(xtxinit);
        //(xT*x)^(-1)*xT
        xtxinvxtinit=erareginit.mult(xtxinvinit,xtinit);
        //(xT*x)^(-1)*xT*y
        betareginit=erareginit.mult(xtxinvxtinit,yreginit);
        return betareginit;
    }    
//------------------------------------------------------------------------------     
    //define the matrix for independent variables
    private double[][] xregfore;
    //define the beta
    private double[][] betafore;
    //define the # of replicates
    private int rowfore;
    //define the # of independent variables
    private int colfore;
    //define the forecasting result
    private double[][] resultfore;
    
    //forecasting
    double[][] regfore(double[][] array1,double[][] array2)
    //input: array1, independent variables; array2, beta
    //Column 1: constant (1); Column 2: cooling load at time t-24; Column 3: cooling load at time t-168;
    //Column 4: daily maximum temperature; Column 5: daily minimum temperature; Column 6: temperature at time t;
    //Column 7: cooling load at time t-24; Column 8: cooling load at time t-168; Column 9: relative humidity at time t;
    //Column 10: relative humidity at time t-24; Column 11: relative humidity at time t-168;
    //Column 12: temperature square at time t; Column 13: temperature square at time t-24;
    //Column 14: temperature square at time t-168; Column 15: product of temperature and relative humidity at time t;
    //Column 16: product of temperature and relative humidity at time t-24;
    //Column 17: product of temperature and relative humidity at time t-168;        
            
    //output: forecasted value
            
    {
        //function call for matrix operation
        mxop erafore=new mxop();
        //obtain the demension of input matrix
        rowfore=array1.length;
        colfore=array1[0].length;
        //copy the independent variables
        xregfore=new double[rowfore][colfore];
        for (int i=0;i<rowfore;i++)
        {
            System.arraycopy(array1[i], 0, xregfore[i], 0, colfore);
        }
        //copy beta
        betafore=new double[colfore][1];
        for (int i=0;i<colfore;i++)
        {
            System.arraycopy(array2[i], 0, betafore[i], 0, 1);
        }
        resultfore=erafore.mult(xregfore, betafore);
        return resultfore;   
    }
//------------------------------------------------------------------------------   
    //define the maxtric for temporarily saving input
    private double[][][] inputtemp;
    //define the beta for initial regression
    private double[][] betatemp;
    //define the beta for final regression
    private double[][] betaregtemp;
    //define the beta for output
    private double[][] betareg;
    //define the tempprary matrix for regssion
    private double[][] inputregtemp;
    //matrix for temporarily saving independent variablels and dependent variable
    private double[][] xregtemp;
    private double[][] yregtemp;
    //# of replicates of original input
    private int rowreg;
    //# of original input
    private int colreg;
    //matrix for temporarily saving the initial forecast of dependent variable
    private double[][] yforetemp;
    //error of initial forecast
    private double[] difffore;
    //variance of initial forecast
    private double varfore;
    //standard deviation of initial forecast
    private double sdfore;
    //# of replicates without outliers
    private int count;
    //input and dependent variable for final regression
    private double[][] xreg;
    private double[][] yreg;
    //difine the matrix for beta adjustment (base load adjustment)
    private double[][] yregfore;
    //define the adjustment beta
    private double betaadj;
    
    //regression without 3 sigma outliers
    double[][] reg(double[][] array)
    //input: array1
    //Column 1: constant (1); Column 2: cooling load at time t-24; Column 3: cooling load at time t-168;
    //Column 4: daily maximum temperature; Column 5: daily minimum temperature; Column 6: temperature at time t;
    //Column 7: cooling load at time t-24; Column 8: cooling load at time t-168; Column 9: relative humidity at time t;
    //Column 10: relative humidity at time t-24; Column 11: relative humidity at time t-168;
    //Column 12: temperature square at time t; Column 13: temperature square at time t-24;
    //Column 14: temperature square at time t-168; Column 15: product of temperature and relative humidity at time t;
    //Column 16: product of temperature and relative humidity at time t-24;
    //Column 17: product of temperature and relative humidity at time t-168;
    //Column 18: cooling load at time t (dependent variable)                
            
    //output: beta obtainted based on the data without outliers
    {
        //function call for matrix operation
        mr24 proreg=new mr24();
        inputtemp=proreg.reginppre24(array);
        //obtain the # of independent variables
        colreg=inputtemp[0][0].length+3+3; //=17
        //define the matrix of beta
        betareg=new double[48][colreg];
        //calculate beta for the model of each hour
        for (int i=0;i<48;i++)
        {
            betaadj=0;
            //temporarily saving the preprocessed input
            inputregtemp=proreg.regpre24(inputtemp,i);
            //obtain the # of valid replicates
            rowreg=inputregtemp.length;      
            //define the matrix for independent variables, dependent variable, beta for one hour model
            xregtemp=new double[rowreg][colreg];
            yregtemp=new double[rowreg][1];
            betatemp=new double[colreg][1];
            //copy the independent and dependent variables from preprocessed data
            for (int j=0;j<rowreg;j++)
            {
                //copy the independent variables
                System.arraycopy(inputregtemp[j], 0, xregtemp[j], 0, colreg);
                //copy the dependent variable
                yregtemp[j][0]=inputregtemp[j][17];
            }
            //calculate beta through initial regression
            betatemp=proreg.reginit(xregtemp,yregtemp);
                        //calculate the first forecast of dependent variable for outlier removal
            yforetemp=proreg.regfore(xregtemp,betatemp);
            //calculate error and deviation
            varfore=0;
            difffore=new double[rowreg];
            for (int j=0;j<rowreg;j++)
            {
                //absolute value of error
                difffore[j]=yregtemp[j][0]-yforetemp[j][0];
                //variance
                varfore=varfore+difffore[j]*difffore[j];
            }
            //standard deviation
            sdfore=Math.sqrt(varfore);
            //count the # of replicates in which the error is within 3 sigma
            count=0;
            for (int j=0;j<rowreg;j++)
            {
                if (difffore[j]<3*sdfore)
                {
                    count=count+1;
                }
            }
            //save the input and dependent variable in the replicates without 3 sigma outliers
            xreg=new double[count][colreg];
            yreg=new double[count][1];
            for (int j=0;j<rowreg;j++)
            {
                if (difffore[j]<3*sdfore)
                {
                    System.arraycopy(xregtemp[j], 0, xreg[j], 0, colreg);
                    yreg[j][0]=yregtemp[j][0];
                }
            }
            //final regression
            betaregtemp=proreg.reginit(xreg,yreg);
            //copy beta
            for (int j=0;j<colreg;j++)
            {
                betareg[i][j]=betaregtemp[j][0];
            }
            //calculate the second forecast of dependent variable for beta adjustment 
            yregfore=proreg.regfore(xreg,betaregtemp);
            //make adjustment to the beta0 for weekday hours
            if (i<24)
            {   
                //mean shift in past 2 weeks (10 weekdays)
                for (int j=count-10;j<count;j++)
                {
                    betaadj=betaadj+(yregfore[j][0]-yreg[j][0]);
                }
                betareg[i][0]=betareg[i][0]-betaadj/10;
            }
            //make adjustment to the beta0 for weekend hours
            else
            {
                //mean shift in past 2 weeks (4 weekends)
                for (int j=count-4;j<count-1;j++)
                {
                    betaadj=betaadj+(yregfore[j][0]-yreg[j][0]);
                }  
                betareg[i][0]=betareg[i][0]-betaadj/4;
            }
        }
        return betareg;
    }
}





//==============================================================================
//define the class for forecasting
class mf
{
    private int hourfore;
    private int rowfore;
    private int[][] timefore;
    private double[][] xfore;
    private double[] resultfore;
    private int nfore;
    double[] fore(double[][] array1,double[][][] array2)
    {
        rowfore=array1.length;
        nfore=0;
        resultfore=new double[48];
        timefore=new int[48][2];
        xfore=new double[48][17];
        //data sorting for forecasting
        for (int i=rowfore-168;i<rowfore;i++)
        {
            if (array1[i][8]==1)
            {
                if(array1[i][0]>3 & array1[i][0]<11)
                {
                    if(array1[i][1]==1)
                    {
                        hourfore=(int)array1[i][2];
                        timefore[nfore][0]=(int)array1[i][1];
                        timefore[nfore][1]=(int)array1[i][2];
                        xfore[nfore][0]=1;
                        xfore[nfore][1]=array1[i-72][3];
                        xfore[nfore][2]=array1[i-168][3];
                        xfore[nfore][3]=array1[i][4];
                        xfore[nfore][4]=array1[i][5];
                        xfore[nfore][5]=array1[i][6];
                        xfore[nfore][6]=array1[i-72][6];
                        xfore[nfore][7]=array1[i-168][6];
                        xfore[nfore][8]=array1[i][7];
                        xfore[nfore][9]=array1[i-72][7];
                        xfore[nfore][10]=array1[i-168][7];
                        xfore[nfore][11]=array1[i][6]*array1[i][6];
                        xfore[nfore][12]=array1[i-72][6]*array1[i-72][6];
                        xfore[nfore][13]=array1[i-168][6]*array1[i-168][6];
                        xfore[nfore][14]=array1[i][6]*array1[i][7];
                        xfore[nfore][15]=array1[i-72][6]*array1[i-72][7];
                        xfore[nfore][16]=array1[i-168][6]*array1[i-168][7];
                        nfore=nfore+1;
                    }
                    else if(array1[i][1]<6)
                    {
                        hourfore=(int)array1[i][2];
                        timefore[nfore][0]=(int)array1[i][1];
                        timefore[nfore][1]=(int)array1[i][2];
                        xfore[nfore][0]=1;
                        xfore[nfore][1]=array1[i-24][3];
                        xfore[nfore][2]=array1[i-168][3];
                        xfore[nfore][3]=array1[i][4];
                        xfore[nfore][4]=array1[i][5];
                        xfore[nfore][5]=array1[i][6];
                        xfore[nfore][6]=array1[i-24][6];
                        xfore[nfore][7]=array1[i-168][6];
                        xfore[nfore][8]=array1[i][7];
                        xfore[nfore][9]=array1[i-24][7];
                        xfore[nfore][10]=array1[i-168][7];
                        xfore[nfore][11]=array1[i][6]*array1[i][6];
                        xfore[nfore][12]=array1[i-24][6]*array1[i-24][6];
                        xfore[nfore][13]=array1[i-168][6]*array1[i-168][6];
                        xfore[nfore][14]=array1[i][6]*array1[i][7];
                        xfore[nfore][15]=array1[i-24][6]*array1[i-24][7];
                        xfore[nfore][16]=array1[i-168][6]*array1[i-168][7];
                        nfore=nfore+1;
                    }
                    else if(array1[i][1]==6)
                    {
                        hourfore=(int)array1[i][2]+24;
                        timefore[nfore][0]=(int)array1[i][1];
                        timefore[nfore][1]=(int)array1[i][2];
                        xfore[nfore][0]=1;
                        xfore[nfore][1]=array1[i-144][3];
                        xfore[nfore][2]=array1[i-168][3];
                        xfore[nfore][3]=array1[i][4];
                        xfore[nfore][4]=array1[i][5];
                        xfore[nfore][5]=array1[i][6];
                        xfore[nfore][6]=array1[i-144][6];
                        xfore[nfore][7]=array1[i-168][6];
                        xfore[nfore][8]=array1[i][7];
                        xfore[nfore][9]=array1[i-144][7];
                        xfore[nfore][10]=array1[i-168][7];
                        xfore[nfore][11]=array1[i][6]*array1[i][6];
                        xfore[nfore][12]=array1[i-144][6]*array1[i-144][6];
                        xfore[nfore][13]=array1[i-168][6]*array1[i-168][6];
                        xfore[nfore][14]=array1[i][6]*array1[i][7];
                        xfore[nfore][15]=array1[i-144][6]*array1[i-144][7];
                        xfore[nfore][16]=array1[i-168][6]*array1[i-168][7];
                        nfore=nfore+1;
                    }   
                    else
                    {
                        hourfore=(int)array1[i][2]+24;
                        timefore[nfore][0]=(int)array1[i][1];
                        timefore[nfore][1]=(int)array1[i][2];
                        xfore[nfore][0]=1;
                        xfore[nfore][1]=array1[i-24][3];
                        xfore[nfore][2]=array1[i-168][3];
                        xfore[nfore][3]=array1[i][4];
                        xfore[nfore][4]=array1[i][5];
                        xfore[nfore][5]=array1[i][6];
                        xfore[nfore][6]=array1[i-24][6];
                        xfore[nfore][7]=array1[i-168][6];
                        xfore[nfore][8]=array1[i][7];
                        xfore[nfore][9]=array1[i-24][7];
                        xfore[nfore][10]=array1[i-168][7];
                        xfore[nfore][11]=array1[i][6]*array1[i][6];
                        xfore[nfore][12]=array1[i-24][6]*array1[i-24][6];
                        xfore[nfore][13]=array1[i-168][6]*array1[i-168][6];
                        xfore[nfore][14]=array1[i][6]*array1[i][7];
                        xfore[nfore][15]=array1[i-24][6]*array1[i-24][7];
                        xfore[nfore][16]=array1[i-168][6]*array1[i-168][7];
                        nfore=nfore+1;
                    }
                }
            }
        }
        for (int i=0;i<48;i++)
        {
            if (i<24)
            {
                if (timefore[i][0]<6)
                {
                    hourfore=timefore[i][1];
                } 
                else
                {
                    hourfore=timefore[i][1]+24;
                }
                for (int j=0;j<17;j++)
                {
                    resultfore[i]=resultfore[i]+xfore[i][j]*array2[hourfore][j][0];
                }
                xfore[i+24][1]=resultfore[i];
            }
            else
            {
                if (timefore[i][0]<6)
                {
                    hourfore=timefore[i][1];
                } 
                else
                {
                    hourfore=timefore[i][1]+24;
                }
                for (int j=0;j<17;j++)
                {
                    resultfore[i]=resultfore[i]+xfore[i][j]*array2[hourfore][j][0];
                }
            }
        }
        return resultfore;
    }
}


//==============================================================================
//define the class for matrix calculation operations
class mxop 
{
 //----------------------------------------------------------------------------- 
    //define the variables of column and row for the matrix transpose
    private int coltrans;
    private int rowtrans;
    //define the result of maxtrix transpose
    private double[][] resulttrans;
    
    //matrix transpose
    double[][] trans(double[][] array)
    //input: array, objective matrix
            
    //output: transposed matrix
            
    {
        //obtain the demension of the objective matrix
        rowtrans=array.length;
        coltrans=array[0].length;
        //define the transposed matrix
        resulttrans=new double[coltrans][rowtrans];
        //transpose the matrix
        for (int i=0;i<coltrans;i++)
        {
            for (int j=0;j<rowtrans;j++)
            {
                resulttrans[i][j]=array[j][i];
            }
        }
        return resulttrans;
    }
//------------------------------------------------------------------------------    
    //define the variables of column and row for the matrix multuply
    //column and row for multiplicand
    private int colmult1; 
    private int colmult2;
    //column and row for multiplicator
    private int rowmult1;
    private int rowmult2;
    //result of matrix mupltiply
    private double[][] resultmult;
    
    //matrix multiply
    double[][] mult(double[][] array1,double[][] array2)
    //input: array1, multiplicand matrix; array2, multiplicator matrix
    
    //output: product of two objective matrixes, array1*array2
    {
       //obtain the demension of the two matrixes
       rowmult1=array1.length;
       rowmult2=array2.length;
       colmult1=array1[0].length;
       colmult2=array2[0].length;
       //define the matrix of production
       resultmult=new double[rowmult1][colmult2];
       //calculate each cell of the matrix
       for (int i=0;i<rowmult1;i++)
       {
           for (int j=0;j<colmult2;j++)
           {
               for (int k=0;k<colmult1;k++)
               {
                   resultmult[i][j]=resultmult[i][j]+array1[i][k]*array2[k][j];
               }
           }
       }
        return resultmult;
    }    
//------------------------------------------------------------------------------   
    //define the variables of matrix demention for the matrix inversion
    private int ninv;
    //define the matrix for UL decomposition
    private double[][] Uinv;
    private double[][] Linv;
    private double[] y1inv;
    private double[] y2inv;
    //define the result of matrix inversion
    private double[][] resultinv;
    
    //matrix inversion using UL decomposition   
    double[][] inv(double[][] array)
    //input: array, objective matrix
            
    //output: inversed matrix
            
    {
        //obtain the demension of the objective matrix
        ninv=array.length;
        //define the inversed matrix
        resultinv=new double[ninv][ninv];
        //define the initial matrix of U for UL decomposition
        Uinv=new double[ninv][ninv];
        for (int i=0;i<ninv;i++)
        {
            System.arraycopy(array[i], 0, Uinv[i], 0, ninv);
        }
        //calculate matrix of L for UL decomposition
        Linv=new double[ninv][ninv];
        for (int i=0;i<ninv;i++)
        {
            for (int j=i;j<ninv;j++)
            {
                if (i==j)
                {
                    Linv[i][j]=1;
                }
                else
                {
                    Linv[j][i]=Uinv[j][i]/Uinv[i][i];
                    Uinv[j][i]=0;
                    for (int k=i+1;k<ninv;k++)
                    {
                        Uinv[j][k]=Uinv[j][k]-Uinv[i][k]*Linv[j][i];
                    }
                }
            }
        }
        //calculate the inversion matrix
        for (int i=0;i<ninv;i++)
        {
            y1inv=new double[ninv];
            y2inv=new double[ninv];
            //solving the equation system of L*temp=y
            y2inv[i]=1;
            for (int j=0;j<ninv;j++)
            {
                y1inv[j]=y2inv[j];
                if (j>0)
                {
                    for (int k=0;k<=j-1;k++)
                    {
                        y1inv[j]=y1inv[j]-Linv[j][k]*y1inv[k];
                    }
                }               
            }
            //solving the equation system of U*ainv[][i]=temp
            for (int j=ninv-1;j>=0;j--)
            {
                resultinv[j][i]=y1inv[j]/Uinv[j][j];
                if (j<ninv-1)
                {
                    for (int k=ninv-1;k>j;k--)
                    {
                        resultinv[j][i]=resultinv[j][i]-resultinv[k][i]*Uinv[j][k]/Uinv[j][j];
                    }
                }
            }
        }
        return resultinv;
    }    
}
